Spinodal decomposition of an ABv model alloy: Patterns at unstable surfaces 
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We develop mean-field kinetic equations for a lattice gas model of a binary alloy with vacancies 
(ABv model) in which diffusion takes place by a vacancy mechanism. These equations are applied to 
the study of phase separation of finite portions of an unstable mixture immersed in a stable vapor. 
Due to a larger mobihty of surface atoms, the most unstable modes of spinodal decomposition 
are localized at the vapor-mixture interface. Simulations show checkerboard-like structures at the 
surface or surface-directed spinodal waves. We determine the growth rates of bulk and surface modes 
by a linear stability analysis and deduce the relation between the parameters of the model and the 
structure and length scale of the surface patterns. The thickness of the surface patterns is related 



C/3 ' to the concentration fluctuations in the initial state. 

I. INTRODUCTION 

C^ , The phase separation of alloys and other unstable mixtures (binary fluids, glasses, polymer blends) has been extensively 
studied for decades pHq]. In a typical experiment, a stable mixture is quenched into an unstable state by a sudden 
r i change in an external control parameter, mostly temperature or pressure. Then, domains of the new equilibrium 
■ O ' phases develop and coarsen, leading to a heterogeneous material. As the mechanical and transport properties of an 
alloy may depend considerably on its microstructure, the understanding of the dynamics of domain formation and 
growth is of great practical importance. In addition, phase separation is an interesting example of a process where 
spontaneous pattern formation occurs during the approach to equilibrium. 

In recent years, attention has been drawn to surface effects on spinodal decomposition, mainly due to interesting 
experimental results on polymer blends [|6|-y . Thin films of mixture are placed between glass plates, or between 
a substrate and vacuum, and quenched. If the interactions with the surface favor one of the components, this 
(<— >) ■ component rapidly segregates to the surface and triggers phase separation in the bulk: one observes surface-directed 
fT^ \ spinodal decomposition, during which concentration waves propagate from the surface into the sample. In some of the 
\^ ' experiments, modulations parallel to the surface are observed, leading to a pattern along the surface with a distinct 
^D ' length scale RH] . The surface effects are in competition with bulk spinodal decomposition, and only a thin layer is 
CTN , affected by the surface, whereas in the interior of the sample the usual bulk structures are found. Interesting questions 
are then what determines the type of structures found at the surface, and to predict typical length scales and growth 
rates of the patterns, as well as the thickness of the surface layer. 

The dynamics of spinodal decomposition near planar substrates has been investigated using continuum models 

C [O-H , Monte Carlo simulations |jla] , and mean- field kinetic equations |l6 L On the other hand, very little is known 

i-rt about free surfaces, which can deform during the decomposition process |17|. We have recently found that surface 

(~| ' modes and spinodal waves at free surfaces can be observed in a simple lattice model p8| , and the present paper is 

O ■ devoted to a detailed study of this phenomenon. 

J-^ ■ The classical theory for the early stage of spinodal decomposition, based on out-of-equilibrium thermodynamics, was 

proposed by Cahn and Hilliard jl| . The Cahn-Hilliard equation is obtained by postulating that the local intcrdiffusion 

currents are proportional to the gradients of the local chemical potential, and supposing that this chemical potential 

. _ can be derived from a free energy functional of Ginzburg-Landau type. The proportionality constant is the atomic 

jrt ' mobility, which is a phenomenological parameter in this theory. For a homogeneous initial state, the Cahn-Hilliard 

■ ■ ' equation may be linearized around the average composition, and analyzed in terms of Fourier modes. Long wavelength 

perturbations grow exponentially with time, whereas short wavelength fluctuations are damped by the gradient terms, 

related to the surface tension. A typical length scale of the resulting domain pattern is then given by the wavelength 

of the fastest growing mode. 

It has been recognized that simplified lattice models are valuable tools to study the influence of microscopic dynamics 
on the process of phase separation in metallic alloys. They provide a convenient conceptual framework with a minimum 
number of parameters. Such lattice gas models assume the existence of a fixed crystal lattice, the sites of which can 
be occupied by the atoms of different species, or by vacancies. A configuration evolves by atomic jumps from site to 



site. Many studies have been carried out on the kinetic Ising model with Kawasaki spin exchange dynamics [|19|-p3[, 
which corresponds to the direct exchange of atoms. In most alloys, however, the dominant mechanism of diffusion 
is a vacancy mechanism 124] , and several models of vacancy- mediated spinodal decomposition have been investigated 
|P5|-P8|. In our model, the rates for atomic jumps follow an Arrhcnius law. 

A standard procedure to investigate the dynamics of such models are Monte Carlo simulations W . The advantage 
of this approach is that the simulations constitute a genuine realization of the stochastic model, and all correlations in 
space and time are preserved. But for the same reason, an analytic understanding is difhcult. Therefore, there have 
been numerous attempts to formulate microscopic kinetic equations which are analytically tractable. These equations 
are derived from the microscopic master equation by an approximation of mean-field type. While the linearized 
equation of Khatchaturyan [ p9| is valid only near equilibrium, recently several authors have proposed fully nonlinear 
kinetic equations for Ising and lattice gas models [50-34|. These equations can be cast in the form of generalized Cahn- 



Hilliard equations, where the atomic mobility depends explicitly on the details of the microscopic dynamics. Such 
equations have been applied successfully to the study of phase-separation dynamics in binary and ternary systems 
p5|-p7t, and to dendritic growth |38). In particular, kinetic equations for models with vacancy dynamics have been 
proposed by Chen and Geng ||3^] and by Puri and Sharma |4^] . We will follow the approach developed by one of the 
present authors p2^ , which allows to relate the microscopic equations in a particularly straightforward manner to the 
equations of out-of equilibrium thermodynamics. 

All the studies of vacancy-mediated phase separation we are aware of start from a homogeneous initial state. We 
study here the evolution of droplets of a dense mixture with few vacancies, immersed in a bath of vacancy-rich "vapor" 
phase. The interesting point in this situation is that we can have interfaces between an unstable mixture and a stable 
vapor which are "neutral" , that is without segregation of one of the components. Such interfaces are impossible in a 
binary model like the Ising model, where any stable phase below the critical temperature favors one of the components. 

In our model, the atoms situated near the surface of the droplets have a larger mobility than bulk atoms. Therefore, 
phase separation is fastest at the surface. We observe different structures at the surface, depending on the model 
parameters and the initial composition of the mixture. In the specially symmetric case of equal interaction energies 
and concentrations for both components, we obtain a regularly modulated surface mode which generates an ordered 
pattern in a surface layer. For asymmetrical interaction energies or compositions, surface-directed spinodal waves 
are observed. This difference can be explained by a competition between surface spinodal decomposition and surface 
segregation. We determine the characteristic length scales and growth rates of bulk and surface modes by a linear 
stability analysis of the mean-field kinetic equations. Our approach allows thus to relate morphological parameters 
to the interaction parameters of the microscopic model. 

The remainder of this article is organized as follows: in Sec. 2, we present the model and the derivation of the 
mean-field kinetic equations. Sec. 3 describes the simulations, which are analyzed in Sec. 4: a linear stability analysis 
is performed in the bulk and at the surface. Sec. 5 presents the discussion of the results. 

II. MODEL AND MEAN FIELD KINETIC EQUATIONS 

We consider a simple cubic lattice in d dimensions, with lattice constant a and a total number of N sites. The sites 
can be occupied by atoms of two species, A and B, or by vacancies v. The occupation numbers at site i, nf, where a 
denotes A, B, or u, are equal to 1 if site i is occupied by species a, and zero otherwise. Double occupancy is forbidden, 
which gives the constraint nf + nf +n\ = 1 Mi. Hence there are two independent variables per site; we choose in the 
following nf and nf . We consider only interactions between nearest neighbor atoms: the energy of a configuration 
C — {nf, nf ,i = 1 . . . N} is given by the Hamiltonian 

H{C) = ~Y. ^AAnfnf + esBnfnfeABinfnf + nfnf)] , (1) 

where the sum is over all nearest neighbor pairs {i, j), and eap are the interaction energies between two atoms occupying 
two nearest neighbor sites. Vacancies do not interact with atoms or other vacancies. Unlike in a binary model, where 
the phase diagram is completely determined by the exchange energy saa + ^bb — "^^ab, here the interaction energies 
are independent. One of them sets the temperature scale, the other two are free parameters. It can be shown that 
Eq. (of) is equivalent to the Hamiltonian of a general ternary system (see appendix) . 

The configuration evolves by jumps of atoms to a neighboring vacant site. This is a common diffusion mechanism in 
metals, and the associated energy barrier is usually much lower than for a direct interchange of atoms. We therefore 
will completely neglect the latter process. As appropriate for activated processes, the hopping rates are assumed to 
follow an Arrhenius law. The underlying physical picture is that the atoms are trapped in potential wells located 
around the sites of the lattice. Atoms spend most of their time near a lattice site, but from time to time they jump 



over the energy barrier between neighboring sites, the necessary energy being provided by the lattice phonons. We 
assume two contributions to the barriers: a constant activation energy Ua and the local binding energy, dH/dnf, 
which depends on the local configuration. The jump rate from site i to j is 



w^^n = I'o exp 






(2) 



Here, T is the temperature (constant throughout the system) and k is Boltzmann's constant. The prefactor Vq is 
related to a vibration frequency of the atom in the well (at least in the absolute Eyring regime) . We can absorb the 
constant Ua in a prefactor which sets the time scale, Wq = Vq exp(— [/^/ZcT). In principle, Wq and Wq are different; 
however, simulation results indicate that qualitative results are unaffected by the ratio Wq/wq as long as it is not too 
far from unity [g5|. Therefore, we will take in the following Wq = Wq ~ wq. The jump rates become 
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Here and in the following, summation over a means a sum over all nearest-neighbor sites. 

The jump rates in Eq. (0) define a stochastic process. To describe its time evolution, one may use the master 
equation for the probability distribution ViC^ t) of finding configuration C at time t: 

^^^ = Y. {V{C',t)W{C'X) - V{C,i)W{C,C')) . (4) 

The transition rates W{C, C) from configuration C to C are equal to w"^ given by Eq. (0) if the two configurations 
differ only by an exchange of a particle and a vacancy, and zero otherwise. Given a solution to the master equation, 
one can formally define a time-dependent average for any operator 0{{nf}), function of the occupation numbers, by 

{0)it)=Y,Oi{nf})P{C,t). (5) 

c 

In particular, we define the time-dependent occupation probabilities of the sites: 

P?{t)^{nf)it). (6) 

These probabilities can also be interpreted as local concentrations. To obtain a kinetic equation for pf(t), one 
differentiates both sides of the above equation with respect to time and uses the master equation. The result is a 
conservation law for the local occupation probability, 

^ = - V i" (7) 

a 

The currents j" through the link i,j are given by 

J- = «(1 - nf - nf)wf_^^ - nf{l - nf ~ nf)w%,) . (8) 

The prefactors of the jump rates assure that the start site is occupied by an a-atom and the target site is empty. 

Up to now, the kinetic equation is equivalent to the complete master equation. To obtain a closed system of equations 
for the occupation probabilities pf , we make a mean- field approximation and replace the occupation numbers in the 
expressions for the currents by their averages pf. Clearly, this approximation is drastic: the resulting equations are 
a set of coupled deterministic differential equations. Fluctuations are suppressed, and the mean-field approximation 
does not take into account correlations. Nevertheless, as in the case of static mean-field approximations, one can 
use this method to get qualitative results on the dynamics of microscopic models. Similar approximations have been 
discussed by several authors |po|-p4| . A possibility to improve systematically the simple mean- field equations is the use 
of the path probability method (PPM) devised by Kikuchi |4|]. The equations of the PPM, however, are considerably 
more complicated, and simulations on bulk spinodal decomposition have shown that the use of the PPM in the pair 
approximation leads to the same qualitative conclusions as the simple mean- field approximation |16| . Therefore, we 
will study in the following only the mean- field case. 

We generalize to our ternary model the method developed for binary systems in Ref. |^2[. The equations for the 
currents can be written as a product of a prefactor 5^-, symmetric with respect to the interchange of i and j, and the 
difference of two local terms Cj and Ci, equivalent to chemical activities: 



JT, = -sr, {cf - ct) , (9) 

with 

Sf^ = ^o(l -vt- Pf )(1 - VJ - Pf) (10) 

and 



^» ~ "i A B ^^P 
^-Pf-Pi 






'i+a 



(11) 



This factorization is not unique (see Ref. p3] for more details). Our particular choice makes it straightforward to 
establish a connection to the phenomenological equations of out-of-equilibrium thermodynamics. We define a local 
chemical potential by 

^i'^^kTlnCr. (12) 

Then, the current becomes 

J- = -Adtj (/i" - Mf ) , (13) 

where the mobility in the link ij is given by 

The explicit expression for the chemical potentials is 



M" = S"— -. (14) 



At" = -eaA^aPf ~ e^B^aPf ~ Z^aApf ~ ZEaBpf + kT In j tt , (15) 

1 - pf - Pi 

where Aagi = X]a(5*+a ^ 9i) for any site-dependent quantity g is the discrete Laplacian, and z is the coordination 
number of the lattice. This expression can also be obtained as the derivative of a free energy function F with respect 
to the local occupation, 

fif = dF/dp1. (16) 

The free energy function is a discrete analog of a Ginzburg-Landau functional: 

F = Y.f{pt.pf,T) + Y: {^ipA.-pff + ^(Pf+. -pf )(pf+. -pf )^(pf+a -pf )') (17) 

i a 

with a local free energy density 

J-/ A B rr\ Z€AA A^ A B ZEbB b^ 

/(P ,P ,T) = —p -zeABP P ^P 

+ kT [p^ \np^+p^ \np^ + (1 - p^ - p^) ln(l - p"^ - p^)] . (18) 

This free energy could have been obtained by a simple static mean field approximation of the Hamiltonian. Further- 
more, for a closed system (no currents crossing the boundaries), F can only decrease. This can be seen explicitly by 
taking its time derivative, using ( |l3[) and noticing that the mobility is always positive: 

i,a ' {i,j),a 

Therefore, the dynamics leads to a state which minimizes the static mean-field free energy. This final state may be the 
ground state (global minimum) or a metastable state (local minimum); we cannot describe nucleation events, unless 
we explicitly introduce fluctuations, for example by adding Langevin noise to the deterministic equations. 

Stated in terms of chemical potentials and mobilities, our kinetic equations have the form of generalized Cahn- 
Hilliard equations. In contrast to the phenomenological equations, the mobilities depend on the local configuration 
and are related to the details of the microscopic jump processes, albeit in an approximate manner. This allows in 
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FIG. 1. Cut through the phase diagram at constant temperature kT — O.GeAA- The interaction energies are eaa = ^bb and 
eAB ~ ^aa/'2- White: one phase, hght gray: two phase coexistence, dark gray: three phase coexistence. 



principle an application of these equations to situations far from equilibrium, and to situations where the mobility 
depends strongly of the local concentration. For the vacancy dynamics considered here, the concentrations of A and 
B are independent variables, but their evolution is coupled by the vacancy field. The mobilities are higher in regions 



with high vacancy concentration. To see this, consider the mobilities for a homogeneous system. For pf 
pf ~* P^j Eq. (p^) becomes 
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(20) 



This expression has a simple interpretation: the prefactor is the mean-field probability of finding an a-atom and a 
vacancy on neighboring sites, and wq times the exponential term is the mean jump rate. 

It is worth noticing that in the mean-field approximation there are no "off-diagonal terms" in the mobility matrix: 
in general, one should obtain an expression for the currents of the form j" = — M'^'^V/i"^ — M"-^V/i^. The reason 
for these terms missing is probably the suppression of all correlations in Eq. (0) by the mean-field approximation. In 
the present context, this deficiency seems to be of little importance. It should be noted that even if the off-diagonal 
mobilities are zero, the same is not true for the diffusion coefficients, because the chemical potentials involve the 
concentrations of both species. 

Prom the free energy density, we can obtain the phase diagram. We will consider only completely attractive models 
here, where no order-disorder transitions occur [caa > 0, cbb > 0, caa + ^bb — '^^ab > 0). The conditions for phase 
coexistence are that the chemical potentials and the grand potential, fi = / — ^^p^ — ji^p^ be equal in the two (or 
three) phases. This is equivalent to a "common tangent plane" construction: the free energy density, function of two 
variables, defines a surface in the space {p^,p^ , /). The (homogeneous) chemical potentials as functions of p'^ and p^ 
define the orientation of planes tangent to this surface. The condition of equal grand potential implies that the points 
representing two (or more) phases in equilibrium must lie in the same plane. Thus we can obtain the phase diagram 
by constructing all the planes that are tangent to the free energy surface in at least two points. Various structures 
of phase diagrams can be obtained p2|. For attractive interactions, quite generally the free energy surface has three 
minima for sufficiently low temperatures. Then, there exists exactly one plane which is tangent in three points: we 
have three-phase coexistence. Besides, there are families of double tangent planes which give the coexistence lines for 
two-phase coexistence. For saa = ^bb, out model is equivalent to the Blume-Emery-Griffiths model ^^. We will 
focus in this paper on the specific example cab = ^aaI^- As we then have tAA = ^bb = ^aa + ess — 2eAB, the phase 
diagram is completely symmetric with respect to the interchange of any two components and may be calculated using 
the analogy with the three state Potts model (see appendix). The phase diagram for kT/eAA = 0.6 and z — 4 (two 
dimensions) is shown in Fig. w. There is a large region where an A-rich, a B-rich, and a vacancy-rich "vapor" phase 
coexist. In the regions of two-phase coexistence, the concentration of the third component is very low. 

Prom formula (EG) we can immediately deduce that the mobility in the dilute "vapor" -phase, where p"^ and p^ 
are small, is much larger than in the two dense phases. As we have caa = ^bb > ^ab, the diffusion of the minority 
component is always faster than that of the majority component in the dense phases. 




FIG. 2. Snapshot pictures of the decomposition of a droplet on a 128 x 128 lattice with periodic boundary conditions. 



Parameters are: €aa = ^bb = 1, ^ab = 
blue, and green are proportional to p^, 
and 149200 10,7 \ 



0.5, kT = 0.5 



Psol 



Psol 



0.46464, and pC 






0.03536. The intensities of red. 



p", and p", respectively. The pictures are taken at t = 14, 24400, 49400, 74400, 99400, 



III. SIMULATIONS 



The part of the phase diagram which is most appropriate for the description of an aUoy is the region of AB-coexistence 
with low vacancy concentration. All studies of lattice gas dynamics with vacancies we are aware of are limited to 
this area. We will investigate the behavior of finite "droplets" of such a material immersed in a stable "vapor" , a 
very natural situation which can arise for instance when droplets of a liquid mixture in coexistence with its vapor 
are rapidly quenched into an unstable state. Evidently, our model is not adapted to describe diffusion processes in a 
vapor, where diffusion does not take place via activated jumps to nearest neighbor sites. But the important feature 
of this "vapor" phase is that diffusion is much faster than inside the "solid" . Moreover, the atomic mobility decreases 
continuously across the vapor-mixture interface, and hence at the surface the diffusion is faster than inside the bulk. 
As we shall see, this induces fast surface modes. 

We integrated the mean-field kinetic equations by an explicit Euler scheme. The time step is limited by the 
numerical stability of the algorithm. For our inhomogeneous system, the "most dangerous" regions of the simulation 
domain are those where the diffusion is fastest, thus the vapor phase. For the diffusion equation with a diffusivity 
-Dvap, the maximum time step is Dva,p/2da^ in d dimensions. For the temperatures and vapor compositions we used, 
Dvap/ol^wo is slightly less than unity. We integrated with a maximum time step of l/4wo in 2D and I/Gwq in 3D 
without encountering numerical instabilities. When we started with step functions as initial conditions, the time step 
had to be chosen much smaller at the beginning and was then slowly increased to the maximum value. 

As initial state, we chose "droplets" (i. e. circular domains) or slabs of a mixture with few vacancies (typically 
some percents) and concentrations p^j and p^j, immersed in a vapor of concentrations p^^p and p^^p. To trigger the 
phase separation, we added small fluctuations to the initial state. To assure mass conservation, pairs of neighboring 
sites and a component (A or B) were randomly chosen, and the concentrations at the two sites were shifted by +Aor 
and —Aqt, respectively, where r is a random number uniformly distributed between —1 and 1. The noise amplitude 
Aq ranged between 10~^ and 10^^. For 2D simulations, we used a lattice of size 128 x 128 with periodic boundary 
conditions. All our simulations were carried out on workstations and took from 2 to 20 hours of CPU time. We also 
simulated two 3D-samples on a 32 x 32 x 64 lattice; these took up to 100 hours CPU time. 

Let us first discuss critical quenches {p^oi — pfoi)- Snapshot pictures from the time evolution of a droplet are 
shown in Fig. 0. Phase separation starts at the surface: a fairly regular modulation appears along the mixture- vapor 
interface. This surface mode triggers phase separation in adjacent regions and propagates into the interior of the 
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FIG. 3. Plots of p„ (a) and p^ (b) along a line normal to a flat interface. From the initial step profile, concentrations rapidly 
relax to a smooth profile. Then, on a longer time scale, phase separation of A and B occurs, starting from the surface. Times 
in units of Wq ^ . 
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FIG. 4. Vacancy distribution during the simulation of Fig. ti. Dark areas have low vacancy content. The contrast inside the 
droplet has been enhanced. Vacancies accumulate in the domain boundaries. The pictures correspond to snapshots number 
3,4, and 6 in Fig. | (f = 49400, 74400, and 149200 w^^). 



sample with a constant velocity, leaving behind a checkerboard-like ordered structure. The domains often coalesce to 
form stripes. 

To see more in detail what happens at the interface, we plot in Fig. several snapshots of the concentration 
profiles along a line which is normal to the surface at a randomly chosen point. The initial step profile quickly relaxes 
to a smooth shape and stays nearly stationary, until on a slower time scale B is enriched and A depleted; at other 
interface points, the opposite happens. An oscillatory concentration profile develops, with an amplitude which has 
its maximum in the interface and decays into the solid. We will show below that the envelope of this oscillation 
becomes a decaying exponential away from the surface. The fact that the perturbation decays exponentially with 
the distance from the surface, but grows exponentially with time, explains the constant propagation velocity of the 
decomposition front. Note that the concentrations in the vapor vary only very slightly: the vapor is a stable phase, 
and hence perturbations decay. Because of the fast diffusion in the vapor, the surroundings of the droplet act as a 
particle reservoir. 

As the surface mode develops, the mixture-vapor interface deforms, and fingers of vapor start to grow into the 
interior of the droplet. This is the result of a MuUins-Sekerka instability E3] with respect to the vacancies. To clarify 
this point, we show in Fig. ^ snapshots of the vacancy concentration during the decomposition process. At the 
beginning, the vacancy concentration stays constant and equal to the initial value. Once the decomposition process 
reaches its nonlinear stage, however, vacancies are expelled from the domains of the new equilibrium phases. These 
excess vacancies have to diffuse to the surface of the droplet. A finger of vapor protruding into the mixture enhances 
the concentration gradients around its tip, and hence grows faster than a flat portion of the surface. In addition, the 
diffusion is faster in the initial mixture than in the phase-separated domains. The growth of the fingers stops when 
a layer of decomposed material has formed around their entire contour. The diffusion then takes place mainly along 
the domain boundaries, where the vacancies are enriched. The fingers are smoothed out by the subsequent coarsening 
process. 
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FIG. 5. Plots of S^{ky,x, t) (a) for the first solid layer and (b) in the middle of the slab. The data are averaged over 10 runs 
with different initial fluctuations (noise amplitude Ao = 0.5 x 10~^). 



The propagation of the surface mode stops when the bulk modes enter the nonlinear regime. The structures in the 
interior of the sample are the usual bicontinuous patterns of bulk spinodal decomposition at equal volume fractions. 
Both surface and bulk structures coarsen by the evaporation-condensation mechanism. This process is faster at the 
surface because of the rapid diffusion through the vapor. At the exterior surface of the droplet, there exist trijunction 
points where the three phases arc in contact. The angles between the interfaces at these points are fixed by the local 
equilibrium between the three surface tensions of Av, Bv, and AB-interfaces. For our symmetric choice of interaction 
energies, all angles are 120° in local equilibrium. 

The initial values for the concentrations in our example represent a special choice: the chemical potentials of the 
two species and the grand potential have the same value in the mixture and in the vapor. There exists exactly one 
set of concentrations satisfying these conditions at a given temperature. In this special case, there is no net mass 
flux between vapor and mixture, and the interface is at rest. This choice was mainly made to simplify the stability 
calculations to be presented below. The existence of a surface mode, however, is not limited to this special case. 
For other initial conditions (but still p^j = p^j and p^p — p^p), a surface mode develops while the mixture- vapor 
interface slowly moves. 

To obtain more quantitative information about the phase separation process, we repeated our simulation with the 
same parameters, but this time in a stripe geometry: a slab of mixture along the y-direction is immersed in the vapor; 
we thus have two straight interfaces normal to the x axis. A convenient quantity for the analysis of phase separation 
processes is the dynamical structure factor. Of particular interest in our case is the difference between bulk and 
surface behavior. Therefore, we define a one-dimensional structure factor along the interface: 
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where Xj and yj are the coordinates of lattice site j. The sum goes over a lattice plane (in 2D: a line) at a fixed 
x-coordinate, and ky is parallel to the surface. Figure ra shows plots of this quantity for different times. We observe the 
characteristic Cahn-Hilliard behavior: in the beginning, linear superposition is valid, and the structure factor for each 
mode grows exponentially with a growth rate depending on ky. In particular, all the structure factor curves intersect 
in one point, corresponding to the marginally stable mode. At later times, the nonlinear terms in the equations of 
motion couple the different modes, and the shape of the structure factor curve changes, as can be noted at the last 
time for the surface layer. The two sets of curves have very different amplitudes, and the maximum is located at 
aky ~ 0.3 at the surface, and at aky « 1 in the bulk. These differences between bulk and surface behavior will be 
addressed in Sec. 4. 

Very different structures appear when the concentration of the mixture is sufficiently off-critical. Fig. g shows the 
evolution of a slab of an AB mixture with a concentration ratio 60:40. The initial concentrations were again chosen 
to give equal chemical potentials and grand potential in the two bulk phases; note that now /i"* ^ ^^ . The minority 
component rapidly segregates at the surface, triggering a "spinodal wave" normal to the surface. A surface mode 
with modulations along the interface is still present and leads to a destabilization of the first layer of the minority 




FIG. 6. Snapshot pictures of a mixture with PsoL_= 0.56007, p^oi = 0.37338 in a vapor of Pvap = 0.04507, Pv 
The other parameters and the colors are as in Fig. § t = 10600, 60700, 123200, 248300, 373400, and 623500 it)(7 



0.03019. 



component, which disintegrates into regularly spaced droplets. These droplets coarsen rapidly. This time, there is no 
MuUins-Sekerka instability with respect to the vacancies, because the formation of the first decomposed layer rapidly 
blocks the exchange of vacancies between the interior of the sample and the vapor. 

Inside the sample, the surface-induced wave travels until the bulk modes reach their nonlinear regime. There is 
a competition of droplet and stripe patterns during the subsequent coarsening process. The droplets in the interior 
of the sample coarsen more rapidly than the stripes at the surface, not surprisingly as the driving force for the 
evaporation-condensation mechanism of coarsening is the curvature of the domain walls. Ultimately, the stripes start 
to break up and are "infected" by the droplet pattern. 

A simulation for a "droplet" geometry is shown in Fig. 0. The obviously symmetric configuration of the outer 
domains of B-rich phase is due to the anisotropy of the surface tensions introduced by the lattice pSl . It is interesting 
to note that this effect is not immediately visible in the critical quenches. The symmetric configuration, however, is 
only transient: on the last snapshot picture, the smallest of the four outer B domains is about to evaporate. Note 
also the symmetry in the first and second ring of inner droplets; this pattern is later destroyed by the coarsening of 
the bulk structures. 

Figure shows the evolution of the interface profiles. In contrast to the critical quench, there is no smooth, nearly 
stationary state at intermediate times. The profile immediately starts to show the onset of oscillations. Also, the 



evolution is faster: at t = 12000z/;q , when the separation of A and B is still tiny in Fig. 



it is already well 



pronounced in the off-critical case. 

The difference between the two evolutions can be qualitatively understood by a look on the free energy surface, 
plotted in Fig. ^j. In this figure, the points a, b, and c denote the initial compositions of the critical and off-critical 
mixtures and the vapor, respectively. In the symmetric case, we can connect the points a and b along a symmetry 
axis of the free energy density, with zero slope along the A-B-direction. This means that an interface which is 
situated completely on this line is stationary with respect to a separation of A and B. For the case of an off-critical 
concentration, we cannot find any such line going from b to c, and an unstable stationary interface does not exist. At 
the surface of the mixture, chemical potential gradients will always lead to a segregation of one of the components to 
the surface. Which component is attracted to the surface depends on the choice of concentrations in the vapor phase. 
Let us mention that this argument is not entirely complete, because it considers only the free energy density, whereas 
the complete free energy also contains the discrete gradient terms. In our simulations, however, we never observed 
any unstable stationary interface configuration for off-critical compositions (or asymmetric interaction energies). 

The transition from the checkerboard structures to the stripes is gradual. For slightly asymmetric concentrations, 
the segregation to the surface is slow, and the surface mode has enough time to grow. We observed checkerboard 









FIG. 7. Same parameters as the preceding figure, but in circular geometry. The symmetries in the picture are a consequence 
of the lattice structure. Times: 10600, 35600, 123200, 248300, 373400, and 623500 Wq^ 
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FIG. 8. Same as Figure H, but for an off-critical quench. Note the faster evolution. Times in units of Wq . 
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FIG. 9. Free energy surface, represented by its isoenergy lines, for kT/eAA = 0.5. The three minima near the corners are 
separated by a "hill" in the center. 

structures up to a concentration ratio of approximately 54:46. For this coniposition, checkerboard structures and 
stripes appear simuhaneously on different portions of the surface. Similar findings are also valid when we vary 
the interaction parameters in our model: for cbb t^ ^aa, we have always observed surface-directed spinodal waves; 
however, a surface mode should appear in this case if for some asymmetric compositions the surface segregation 
becomes slow. 

These findings are consistent with calculations for mixtures near flat substrates using continuous equations of Cahn- 
Hilliard type: spinodal waves occur when one component is attracted to the substrate [Ijjl^, whereas surface modes 
have been found in the case of a substrate which prefers neither of the components of the mixture p4| . Our simple 
model shows that this general behavior is also valid for free surfaces. 

In 3D samples we also find fast surface modes, but because the vapor-mixture interface is now two-dimensional, 
the patterns occurring at the surface for critical quenches are those of 2D bulk spinodal decomposition (Fig. |l^) . A 
bicontinuous pattern forms at the surface and propagates into the sample, replicating itself in an oscillatory manner, 
that is we always find an oscillating concentration profile when we look normal to the surface. In the off-critical case, 
a spinodal wave occurs as in two dimensions (Fig. |l l| ) . 

In all the simulations presented so far, the domains of the A- and B-rich phases stick together, because the surface 
tension of an Av- or Bv-interface is more than the half of the one of an AB-interface. This changes if we lower the 
interaction energy eab'- for tAB = 0, the droplet "explodes": thin layers of vapor penetrate into the interior of the 
droplet along the forming AB-interfaces, and the domains of A and B are slowly drifting apart (Fig. |lj). This reminds 
of the decomposition of a binary mixture in the presence of a surfactant. The presence of sharp corners and facets in 
the domain shapes of the last snapshot indicates that the surface tension anisotropy is quite large. 

IV. LINEAR STABILITY ANALYSIS 

A. Bulk 

We will now study more in detail the checkerboard structures. To this end, we must calculate the growth rates of 
bulk and surface modes. We start with the bulk modes. A homogeneous system of overall composition p"^ and p^ is 
perturbed by small fluctuations of the occupation probabilities: 

pt{t)=p^ + 6f{t) ia = AB), (22) 
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FIG. 10. 3D sample on a 32 x 32 x 64 lattice at critical composition (initial concentrations as in Fig. H, kT/eAA = 0.75). 
The outer white spheres indicate the limits of the simulation box. Periodic boundary conditions are applied in all directions. 
The size of each sphere inside the box indicates the total concentration {p +p^), the color its composition (red: A, blue: B). 
Regions of vapor are situated above and below the mixture film; as the concentration of atoms is very low in these regions, the 
corresponding spheres are tiny. This snapshot was taken at i = 20400u)^^. 




FIG. 11. 3D sample at off-critical composition. Initial concentrations are as in Fig. pi kT/eAA ~ 0.75, and t — 20400u;o 
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with (5f ^ 1. We linearize the chemical potentials around the average concentrations. Introducing a vector notation 



with respect to the two species of particles, we obtain starting from Eqs. (15) 



:$)<^)M^)Mr 



B ,-\ r.B ,-^ aVb +S ,fe • (23) 



Here, ^" are the unperturbed values of the chemical potentials, the matrix E of the interaction energies is 

^^(^AA eAB \ ^24^ 

V <^AB esB J 

and S is the matrix of the second derivatives of the free energy density, taken at p and p^ : 



Saa Sab \ •+!, c _ ^ -^ 
Sab Sbb ) ^^*^ '^"'^ " d^^^ 



(25) 



The variations in the chemical potentials create currents. To obtain these currents to order one in Sf, we may use 
the unperturbed values of the mobilities, 



The equations of motion become: 



M'; = MS,Jp-',p'') = M'^ (26) 



dt \Sf J \ M^Aa^i 



As a homogeneous system is translation invariant with respect to the lattice vectors, solutions of the linearized 
equations are of the form: 

where Xj is the position vector of site j in real space, and k = (kx, ky) is the wave vector of the perturbation. The 
growth rate uj{k) and the coefRcients 5°' have to be determined by solving the eigenvalue problem 

-(^1 f 'b ) = A, { ^; ,-0, ) (S - A,E) (f,). (29) 

Here, the terms 



5" / MO A/" / ^ k-J I s 



Aj: = -4sin2(fc^a/2) -Asm^{kya/2) (30) 

arise from the discrete Laplacians. Eq. ( p9| ) is quadratic in cu and thus gives a stability spectrum with two branches. 
Each of these branches can be stable (w is negative for all wave vectors) or unstable. In the latter case, positive growth 
rates occur for small values of |fc|. The number of unstable branches is equal to the number of negative eigenvalues of 
the matrix S. This can be easily seen by taking the limit \k\ ^ in Eq. (|29|). The terms proportional to Ai can be 
neglected, and we have 



c^{k)(^lt)^-\k\'s(^lty (31) 

The matrix S is related to the curvature of the free energy surface. If both eigenvalues are positive, the surface is 
locally convex, and the homogeneous state is stable. For two eigenvalues of different sign, the surface has locally the 
structure of a saddle point, and we have partial instability: only fluctuations in the concave direction in concentration 
space are amplified. Finally, for two negative eigenvalues, the free energy surface is concave and all perturbations 
grow. The frontiers between these regions of different stability behavior are given by the spinodal surfaces, defined 
by the condition 

det S = 0. (32) 
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This generalizes the concept of a spinodal curve to our three-component system. 

For the special case of symmetric interaction energies and equal average concentration, Eq. (|29| ) can be considerably 
simplified because then we have M^ — M^ — M. In addition, the matrix L defined by 

L = ^,,(S - ^gE) (33) 

is symmetric. We obtain immediately 

uj±{k)=M(LAA±LAB), (34) 



and the associated eigenvectors are 



SB ]-\±i 



B I - I 11 I • (35) 



This last result shows that of the two branches of the dispersion relation, the first describes the separation into a 
dense AB mixture and a dilute vapor, whereas the second gives a separation between A and B, leaving the local 
vacancy concentration unchanged. Fig. |l^ shows a cut through the spinodal surfaces along the axis p^ = P^ ■ We also 
indicate the regions of different stability behaviors, and show the typical shapes of uj+{k) and uj-{k) in these regions. 
The interior of our samples is in the region where only the mode separating A and B is unstable: the local vacancy 
concentration stays unchanged in the linear stage. Other modes of decomposition in a homogeneous ternary system 
have been studied by Chen | p5[ |. 

The anisotropy due to the lattice structure enters in the above formulae by the factor Aj^. For the long wavelength 

perturbations considered here {a\k\ < 1), the relative variations of Aj^ with orientation are of the order of a percent. 
We will therefore neglect this dependence and use a wave vector along one of the lattice directions for comparisons to 
numerical results. For the example of the simulation shown in Fig. g, i.e. with kT /caa — 0.5 and p"^ = p^ = 0.46464, 
we found a maximum bulk growth rate of ojb — 1.06 x 10~'*u>o at a wave number kb ~ 1.003/a. 

B. Surface 

We must now analyze the stability spectrum at the surface. This is more difficult than in the bulk because the initial 
state is now heterogeneous. To simplify the problem, we will treat the case of a flat interface which is normal to one of 
the lattice directions, say x. Then, in the initial state all lattice sites in a layer at a given x coordinate have the same 
concentrations. In what follows, we will replace the site indices "i" used so far by a pair of indices "n, j", where n 
numbers the layer, and j numbers the sites in the y direction. We use the same values of the concentrations in the bulk 
phases as for our simulations. As the chemical potentials are equal in the two phases, the initial interface state can 
be obtained by fixing the chemical potentials to their appropriate values and numerically solving the one-dimensional 
version of the finite difference equations Eq. ( [l5| ) for the concentrations in the nth layer p" . 

Whereas the translation invariance is broken along the x-axis, it is preserved along y, and hence we can still use a 
Fourier representation. We write the perturbed state as 

PnAt)^P?.+Sl,{t), (36) 

SnA^) = ^n exp (iky J a + Loiky)t) , (37) 

with two unknown coefficients (5„ and S^ per layer. We start with the linearization of the chemical potentials. The 
discrete Laplacians appearing in Eq. ([l5| ) give contributions of the form 

A A" — A" -I- A" ^ A" 4- A" — zLA" 

= [-5^1 + ^,'^+1 + {Af: - 2) S"^] exp {ikyja + u;{ky)t) . (38) 

Denoting by p" the (constant) values of the unperturbed chemical potential, and by S(n) the matrix of the second 
derivatives, taken at p„ and p^ , we find 



f^nj - P"" 



S^- \ 




^ (5^ 1 ■ \ 




^ S\, ■ 


"J 1 


E 


n-1,7 1 


E 









(S(n) - {Aj: 2)E) J^^ - E Jif''' - E J}^'^^ . (39) 
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FIG. 14. Maximum growth rate iOa and corresponding wave vector k^ for the vapor-mixture interface with initial conditions 
as in Fig. Q as a function of the number A'^ of layers around the interface included in the eigenvalue problem. 

The next step is the hnearization of the equations of motion. As in the homogeneous case, the mobilities can be 
taken in the unperturbed initial state. But because the concentrations vary through the interface, we cannot use the 
expression for the homogeneous system. The appropriate expressions can be obtained from Eq. (]l4) by taking the 
limit /ij -^ [ii with pj ^ pf. 






Wo- 



{l ~ Pi - p^)il ^ Pi 



Pn 



kT 



■ exp 



kT 



(40) 



(the two indices of Mnm are two layer indices) . The linearized equation of motion becomes 



dt "'^ 



M" 



U^: 



nn yr^nj-i-l 



^Mnj + /^nj-l) + -^nn+1 (^n+lj ^ A^nj) + ^^nn-1 (/^n-lj ~ ^n,j) • 



(41) 



Inserting Eqs. (^) and (^, one obtains an eigenvalue problem for w, the eigenvectors given by the set of unknown 
coefficients S". This means that the matrix to diagonalize is of size 2L x 2L, where L is the number of layers in the x 
direction, and that the corresponding dispersion relation has 2L branches. This is a standard linear algebra problem 
and can be solved numerically. 

In the initial stage of the development, the surface mode is localized in a small number of layers around the 
surface. Hence we can simplify the problem by considering a small "solution region" around the interface. We assume 
perturbations outside a narrow region of N layers centered around the interface to be zero. We diagonalized the 
resulting 2N x 2N matrix numerically and obtained the maximum growth rate, w^, and the corresponding wave 
vector kg. As shown in Fig. |lj, these quantities become independent of iV if a sufficient number of layers is included 

11- 

To compare the results of this stability analysis to the simulations, we extracted the stability spectrum of the surface 
from the structure factor data of Fig. |[ For an exponentially growing Fourier mode, we have S{t) oc exp(2a;i), and the 
growth rate can be obtained as lu — \n(S(t2) / S{ti)) /2(t2 — ii). Fig. |l^ shows that the agreement with the theoretical 
prediction is excellent, except for very small wave numbers. This discrepancy can be explained by the fact that for 
large wavelengths, more layers have to be included, as the diffusion field in the vapor will be appreciably modified up 
to a distance from the interface which is comparable to the wavelength. In conclusion, the maximum growth rate and 
the corresponding wave vector can be accurately predicted by a calculation involving only a small number of layers. 

For the example of Fig. y, we find a maximum growth rate of LOg = 4.47 x 10~'^'Wq at a wave number kg = 0.349/a. 
This surface mode grows more than four times faster than the fastest bulk mode, and at a wavelength about three times 
larger than the wavelength of the typical bulk pattern. We verified indeed in our simulations that the "incubation 
time" , i. e. the time after which the growing perturbations become visible in a given visualization, was about four 
times larger in the bulk than at the surface. 

Having verified that our method gives accurate results, we can apply it to study the behavior of bulk and surface 
modes for different temperatures. The results are given in Table |[ All calculation were performed with TV = 16, except 
for the highest temperature where A^ = 32 was used. The bulk growth rate has a maximum around kT/eAA = 0.6 
and decreases as the temperature is lowered, because the activated dynamics lead to small values of the mobility for 
low temperatures. On the other hand, for the highest temperature the system is near a spinodal surface, and the 
driving force for phase separation is small, leading also to a lower growth rate. The surface growth rate follows the 



16 



0.0004 



S 0.0002 



0.0000 




FIG. 15. Comparison of theoretical prediction (line) and simulation results (circles) for the stability spectrum of the va- 
por-mixture interface with the same parameters as in Fig. bl The theoretical curve is calculated with TV = 16, the growth 
rates in the simulation were extracted from the data in Fig. pi The scatter in the simulation data is due to the finite number 
of different runs. 



TABLE I. Characteristic quantities for bulk and surface modes in two dimensions for different temperatures 



kT/eAA 



0.4 



0.5 



0.6 



0.7 



Pva,p iPvap 

Psol jPsol 

LOb/wo X 1000 

uJs/wQ X 1000 

UJs/uJh 

akb 
akg 



0.0141 
0.4859 
0.016 
0.383 
23.9 
1.147 
0.312 
3.68 



0.0354 
0.4646 
0.106 
0.447 

4.21 
1.003 
0.349 

2.87 



0.0724 
0.4276 
0.255 
0.342 

1.34 
0.793 
0.421 

1.88 



0.1412 
0.3588 
0.008 
0.003 

0.36 
0.222 
0.132 

1.68 



TABLE 11. Same as table Nl but in three dimensions. 



kT/eAA 



0.6 



0.75 



0.9 



1.05 



/^vap 5i^vap 
Psol jPsoI 


LJi/lUo X 1000 


w,/«;o X 1000 


Uaj^b 


akb 


aks 


kb/ks 



0.0141 
0.4859 
0.024 
0.976 

40.7 
1.453 
0.414 

3.51 



0.0354 
0.4646 
0.159 
0.773 

4.86 
1.259 
0.410 

3.07 



0.0724 
0.4276 
0.383 
0.528 

1.37 
0.986 
0.505 

1.92 



0.1412 
0.3588 
0.012 
0.007 

0.61 
0.272 
0.204 

1.33 
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same trends. At the highest temperature considered, Ws is lower than cui,, which means that a proper surface mode 
does not exist any more. The ratio Wg/wf, increases monotonically as the temperature is lowered. It should be noted 
that the mean-field approximation becomes increasingly inaccurate for low temperatures. Also, the interfaces become 
very sharp, which leads to strong lattice effects. 

It is also straightforward to repeat these calculations for three dimensions. In fact, the phase diagram stays the 
same if we scale the temperature by a factor of 1.5. On the other hand, the interface shapes and the growth rates 
are different, because the balance between the discrete Laplacians and the local terms is altered in Eqs. (|5|) for the 
chemical potentials. In the stability calculations, the magnitude of the stability matrix S is modified whereas the 
terms arising from the Laplacians stay the same. The final results for the growth rates and wave vectors, shown in 
Table ||, are similar to the two-dimensional case. 

C. Propagation of the surface patterns 

The surface mode propagates into the interior of the sample, enforcing its characteristic wavelength in the direction 
parallel to the surface. This behavior can be understood by considering the solution of the eigenvalue equation away 
from the interface. Even if the surface modes are localized at the vapor-mixture interface, they are eigenmodes of 
the whole system and grow everywhere with the same growth rate ujs. But away from the interface, the initial state 
becomes homogeneous, and we can use Eq. (p4), derived for the bulk, if we allow for a decay of the amplitude in the 
direction normal to the interface by introducing a complex wave number kx- As the growth rate and the y-component 
of the wave vector are already known from the analysis of the surface instability, k^ is the only remaining unknown. 
The relevant dispersion relation is Lu^{k), and we obtain the equation 

4 {eaa -eAB)- A^ {Saa -Sab) + '^=Q- (42) 

This equation can either be solved numerically using the exact expression for Aj: or analytically with the approximation 
Aj: « —k'^ — ky. The latter method leads to a biquadratic equation for k^- In both cases, there are four solutions of 
the form 

kx = ±k' ± ik" (43) 

with k' and k" real and positive. The two solutions with negative imaginary part diverge in the bulk and have to 
be discarded. The other two solutions lead to modes proportional to exp{±ik'x + ujgt — k"x), that is modes with an 
oscillatory concentration profile and an envelope which decays exponentially with the distance from the surface, but 
grows exponentially in time. One can define a propagation velocity of the decomposition front by the phase velocity 
of the envelope. 

We compared the values for the wavelength normal to the surface, A^ — 27r/fc', and the propagation velocity obtained 
by this method to our simulations and found good agreement pq | . 

The competition between surface and bulk modes leads to the interesting consequence that the thickness of the layer 
near the surface where the surface-directed patterns prevail depends on the strength of the initial fluctuations. This 
can be seen as follows. For an exponentially growing mode with initial amplitude So, the time to reach a threshold 
amplitude (5„; where nonlinear couplings between different modes become important is 

T^-ln^. (45) 

Lo do 

Now, a surface mode "propagates" into the sample approximately between the time r^ when it is well developed at 
the surface, and the time Tf, when the bulk modes reach the nonlinear stage. The distance d the front propagates is 
therefore given by 

d^v{n-r,) = -Lf^-l)lJ-f, (46) 

k" \LUb J do 

The thickness of the surface structures depends logarithmically on the strength of the initial noise. This reasoning 
applies both to the symmetrical checkerboard structures and the surface-directed spinodal waves and was confirmed 
by our simulations. 
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V. CONCLUSIONS 

We have developed mean field kinetic equations to describe the dynamics of a lattice gas model of a binary alloy 
with vacancies (ABv model) in which diffusion takes place by the vacancy mechanism only. Despite the simplicity 
of the model, we observe a rich variety of phase separation patterns at free surfaces between an unstable mixture 
and a stable vapor. The most spectacular effect is a fast surface mode. It creates ordered patterns at the surface 
with a length scale which is clearly distinct from the characteristic scale of the bulk patterns. This mode appears in 
a small range of parameters around a symmetric point where the mixture- vapor interface is "neutral" , that is none 
of the components of the mixture segregates to the surface. On the other hand, if such segregation occurs and is 
rapid enough, the surface mode is suppressed, and instead surface-directed spinodal waves are observed which create 
a striped pattern along the surface. Both patterns, once they have formed at the surface, propagate into the sample 
over a distance which is related to the difference in surface and bulk growth rates and the strength of the initial 
fluctuations. 

Our approach starts from a minimal model with very few assumptions. Therefore, it cannot be used to model a 
particular experimental situation. But it allows to identify some basic ingredients necessary for the formation of such 
surface structures in a fairly well-defined setting, and we can draw some conclusions which should be generally valid. 

The existence of the fast surface mode is related to the fact that in our model the surface atoms are more mobile 
than in the bulk, an assumption which seems reasonable for many interfaces between a dense and a dilute phase. 
The characteristic growth rates of the surface modes can be calculated by a linear stability analysis starting from 
the initial interface profile. We have shown that accurate results can be obtained by solving the resulting eigenvalue 
problem in a small domain around the interface. The characteristics of the propagation of the decomposition front 
can be obtained by connecting the solution of the surface problem to the bulk solution. This method could be applied 
as well to continuum models of Ginzburg-Landau type by using an appropriate discretization. 

Which type of surface structure occurs is related to the time scales for segregation and surface phase separation. 
If the segregation is slow, surface spinodal decomposition is dominant, otherwise spinodal waves are observed. As in 
a realistic system, the interactions between different species are always different, spinodal waves are the more generic 
pattern. The surface mode should show up, however, in a small range of initial compositions if the interactions are 
not too different. More work is needed to clarify this point. 

A very interesting point is that the thickness of the layer of surface patterns depends on the initial fluctuation 
strength. This is not the case in bulk spinodal decomposition, where a rescaling of the initial fluctuations amounts 
simply to a shift in time. The arguments leading to Eq. (^) are fairly general and should apply to all systems where 
a fast surface mode is in competition with a slow bulk mode. In rapid quench experiments, the initial fluctuation 
spectrum is mainly determined by the temperature before the quench. Therefore, the thickness of the surface layer 
may depend on the initial temperature. In view of the logarithmic dependence of d on the noise strength, this effect 
might be difficult to observe; however, if the initial state of the system is close to a critical point, the fluctuation 
amplitude is a rapidly varying function of temperature. 

Surface effects in spinodal decomposition have been studied recently in polymers [|6H11|. Evidently, our equations 
cannot be applied to polymers. But from our results it is simple to construct a Ginzburg-Landau theory in which 
surface modes occur: it is sufficient to introduce a mobility function which explicitly depends on the density and 
which enhances diffusion on the surface. 

It would be interesting to compare our findings to Monte Carlo simulations which naturally contain the fluctuations. 
It has been shown, however, that short-range lattice gas models do not exhibit linearly superposed, exponentially 
growing modes in the early stages of the phase separation. This can be attributed to the fact that the initial noise is 
too strong for a linearization of the equations of motion to be valid ||46|] . On the other hand, for models with longer 
range interactions the Cahn-Hilliard behavior is restored [^, and we would expect surface modes in such models 
even with stochastic dynamics. 

In summary, our approach allows to explore the rich dynamics of the ABv model and to relate the structures that 
form spontaneously during phase separation to the parameters of the microscopic model. Here, we have only explored 
a small part of the possible behaviors in this model, because we have limited ourselves to attractive interactions. 
There are other interesting questions which could be addressed in the framework of this model and using mean fleld 
kinetic equations, for example the influence of the vacancy distribution on the coarsening behavior, or the interplay 
between phase separation, short range ordering and vacancy distribution. 
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FIG. 16. Phase diagram of the completely symmetric ABv-model. We show the equilibrium concentrations in the a (A-rich) 
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APPENDIX: PHASE DIAGRAM 



Let us first show that the Hamiltonian of the ABv-model can be derived from a general ternary Hamiltonian. We 
assume that the sites of the lattice can now be occupied by three different sorts of atoms, A, B, or C. In terms of 



the occupation numbers n^ 



and n^ , the energy of a configuration is given by 



H 






N 



'ap 



tii n- 



EEa^'o^ 



(Al) 



i=l 



The interaction energies in the ternary system are primed to distinguish them from the e's in Eq. ffl). Here, the /i'p 
are external chemical potentials, equivalent in the spin language to "generalized magnetic fields" acting only on one 
possible spin state. Using the constraint of single occupancy, n^ -\- nf -\- nf = 1 Vi, we can eliminate the occupation 
numbers of one species, say C, and we obtain 



N 



^ = - E [^AAnfnf + eBBnfnf + eABinfuf + nfnf)] - ^ {ii^nf + fi^nf) 



(■'J-> 



1=1 



where the effective interaction energies and chemical potentials appearing in this last equation are 



CAA 


= '^AA + '^CC 


- ^^'ac, 


esB 


= f SB + '^CC 


- ^^'bc. 


^AB 


= '^AB + '^CC 


" '^AC ~ '^BC^ 


M^ 


,A iC 

= A*o -Mo - 


' ^^CC + ^<:AC 


A*f 


iB iC 

=Mo -Mo ' 


- Zt(j(j + Z^BQ 



and 



(A2) 



(A3) 
(A4) 

(A5) 
(A6) 
(A7) 



When the total number of particles of each species is conserved, the chemical potential terms in Eq. (A2) are constants. 
Then, we are back to the Hamiltonian Eq. (|l|). 

Let us now consider the three state Potts model, with ternary interaction energies e^^ = e'g^ = ej^^^ = 1 and 
d^n = for a ^ /?. This gives the effective interactions e^A — ^bb — 2, ^ab = 1; by changing the temperature scale, 



20 



we obtain the values caa = ^bb = 1, ^ab = 0.5 used in the present paper. Furthermore, from the symmetry of the 
Potts model it is obvious that three phase coexistence is possible below the critical temperature at zero magnetic 
fields. This is equivalent, in the ABv model, to 

^J.^ = -zeAAP^ - zeABP^ + kT- ^ 77 = -zeAA/2 (A8) 

and an equivalent equation for /i^. These equations have to be solved numerically to obtain the equilibrium concen- 
trations. The result for the A-rich phase is shown in Fig. ^. The coexistence of A-rich, B-rich, and "vapor" -phases 
terminates at a quadruple point: a fourth minimum in the free energy surface, located at the symmetric point 
p^ = p^ = p" =^ 1/3, develops, and at kT — 1/(2 In 2) w 0.721 the four minima are exactly at the same level. Above 
this temperature, there is a narrow temperature range above which the four minima persist, but now the fourth 
minimum is lowest, and we have three different regions of three-phase coexistence (not shown in the figure). Above 
kT =^ 0.75, only the "symmetric" minimum remains. 
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